home *** CD-ROM | disk | FTP | other *** search
/ PC Media 23 / PC MEDIA CD23.iso / share / prog / newmat / newmat.h < prev    next >
Encoding:
C/C++ Source or Header  |  1995-01-17  |  52.7 KB  |  1,489 lines
  1. //$$ newmat.h           definition file for new version of matrix package
  2.  
  3. // Copyright (C) 1991,2,3,4: R B Davies
  4.  
  5. #ifndef MATRIX_LIB
  6. #define MATRIX_LIB 0
  7.  
  8. #include "include.h"
  9.  
  10. #ifdef NO_LONG_NAMES
  11. #define UpperTriangularMatrix UTMatrix
  12. #define LowerTriangularMatrix LTMatrix
  13. #define SymmetricMatrix SMatrix
  14. #define DiagonalMatrix DMatrix
  15. #define BandMatrix BMatrix
  16. #define UpperBandMatrix UBMatrix
  17. #define LowerBandMatrix LBMatrix
  18. #define SymmetricBandMatrix SBMatrix
  19. #define BandLUMatrix BLUMatrix
  20. #endif
  21.  
  22. #ifndef TEMPS_DESTROYED_QUICKLY_R
  23. #define ReturnMatrix ReturnMatrixX
  24. #else
  25. #define ReturnMatrix ReturnMatrixX&
  26. #endif
  27.  
  28. #include "boolean.h"
  29. #include "myexcept.h"
  30.  
  31. /**************************** general utilities ****************************/
  32.  
  33. class GeneralMatrix;
  34. void MatrixErrorNoSpace(void*);                 // no space handler
  35.  
  36. class LogAndSign
  37. // Return from LogDeterminant function
  38. //    - value of the log plus the sign (+, - or 0)
  39. {
  40.    Real log_value;
  41.    int sign;
  42. public:
  43.    LogAndSign() { log_value=0.0; sign=1; }
  44.    LogAndSign(Real);
  45.    void operator*=(Real);
  46.    void ChangeSign() { sign = -sign; }
  47.    Real LogValue() const { return log_value; }
  48.    int Sign() const { return sign; }
  49.    Real Value() const;
  50.    FREE_CHECK(LogAndSign)
  51. };
  52.  
  53. // the following class is for counting the number of times a piece of code
  54. // is executed. It is used for locating any code not executed by test
  55. // routines. Use turbo GREP locate all places this code is called and
  56. // check which ones are not accessed.
  57. // Somewhat implementation dependent as it relies on "cout" still being
  58. // present when ExeCounter objects are destructed.
  59.  
  60. class ExeCounter
  61. {
  62.    int line;                                    // code line number
  63.    int fileid;                                  // file identifier
  64.    long nexe;                                   // number of executions
  65.    static int nreports;                         // number of reports
  66. public:
  67.    ExeCounter(int,int);
  68.    void operator++() { nexe++; }
  69.    ~ExeCounter();                               // prints out reports
  70. };
  71.  
  72.  
  73. /************** Class to show whether to check for loss of data ************/
  74.  
  75. class MatrixConversionCheck : public Janitor
  76. {
  77.    static Boolean DoCheck;
  78. public:
  79.    MatrixConversionCheck() { DoCheck=TRUE; }          // turn check on
  80.    ~MatrixConversionCheck() { DoCheck=FALSE; }        // turn check off
  81.     void CleanUp() { DoCheck=FALSE; }
  82.     static Boolean IsOn() { return DoCheck; }
  83.     static void DataLoss();
  84.     friend class SkipConversionCheck;
  85. };
  86.  
  87. class SkipConversionCheck : public Janitor
  88. // to turn check off in current block
  89. {
  90.     Boolean LastValue;
  91. public:
  92.     SkipConversionCheck() : LastValue(MatrixConversionCheck::DoCheck)
  93.         { MatrixConversionCheck::DoCheck = FALSE; }
  94.     ~SkipConversionCheck() { MatrixConversionCheck::DoCheck = LastValue; }
  95.     void CleanUp() { MatrixConversionCheck::DoCheck = LastValue; }
  96. };
  97.  
  98. /**************************** class MatrixType *****************************/
  99.  
  100. // Is used for finding the type of a matrix resulting from the binary operations
  101. // +, -, * and identifying what conversions are permissible.
  102. // This class must be updated when new matrix types are added.
  103.  
  104. class GeneralMatrix;                            // defined later
  105. class BaseMatrix;                               // defined later
  106. class MatrixInput;                              // defined later
  107.  
  108. class MatrixType
  109. {
  110. public:
  111.    enum Attribute {  Valid     = 1,
  112.                      Diagonal  = 2,             // order of these is important
  113.                      Symmetric = 4,
  114.                      Band      = 8,
  115.                      Lower     = 16,
  116.                      Upper     = 32,
  117.                      LUDeco    = 64 };
  118.  
  119.    enum            { US = 0,
  120.                      UT = Valid + Upper,
  121.                      LT = Valid + Lower,
  122.                      Rt = Valid,
  123.                      Sm = Valid + Symmetric,
  124.                      Dg = Valid + Diagonal + Band + Lower + Upper + Symmetric,
  125.              RV = Valid,     //   don't separate out
  126.              CV = Valid,     //   vectors
  127.              BM = Valid + Band,
  128.              UB = Valid + Band + Upper,
  129.              LB = Valid + Band + Lower,
  130.              SB = Valid + Band + Symmetric,
  131.              Ct = Valid + LUDeco,
  132.              BC = Valid + Band + LUDeco
  133.                    };
  134.  
  135.  
  136.    static nTypes() { return 9; }               // number of different types
  137.                            // exclude Ct, US, BC
  138. public:
  139.    int attribute;
  140. public:
  141.    MatrixType () {}
  142.    MatrixType (int i) : attribute(i) {}
  143.    int operator+() const { return attribute; }
  144.    MatrixType operator+(MatrixType mt) const
  145.       { return MatrixType(attribute & mt.attribute); }
  146.    MatrixType operator*(const MatrixType&) const;
  147.    MatrixType SP(const MatrixType&) const;
  148.    MatrixType operator|(const MatrixType& mt) const
  149.       { return MatrixType(attribute & mt.attribute & Valid); }
  150.    MatrixType operator&(const MatrixType& mt) const
  151.       { return MatrixType(attribute & mt.attribute & Valid); }
  152.    Boolean operator>=(MatrixType mt) const
  153.       { return ( attribute & mt.attribute ) == attribute; }
  154.    Boolean operator==(MatrixType t) const
  155.       { return (attribute == t.attribute); }
  156.    Boolean operator!=(MatrixType t) const
  157.       { return (attribute != t.attribute); }
  158.    Boolean operator!() const { return (attribute & Valid) == 0; }
  159.    MatrixType i() const;                       // type of inverse
  160.    MatrixType t() const;                       // type of transpose
  161.    MatrixType AddEqualEl() const               // Add constant to matrix
  162.       { return MatrixType(attribute & (Valid + Symmetric)); }
  163.    MatrixType MultRHS() const;                 // type for rhs of multiply
  164.    MatrixType sub() const                      // type of submatrix
  165.       { return MatrixType(attribute & Valid); }
  166.    MatrixType ssub() const                     // type of sym submatrix
  167.       { return MatrixType(attribute); }        // not for selection matrix
  168.    GeneralMatrix* New() const;                 // new matrix of given type
  169.    GeneralMatrix* New(int,int,BaseMatrix*) const;
  170.                                                // new matrix of given type
  171.    char* Value() const;                        // to print type
  172.    friend Boolean Rectangular(MatrixType a, MatrixType b, MatrixType c);
  173.    friend Boolean Compare(const MatrixType&, MatrixType&);
  174.                                                // compare and check conv.
  175.    Boolean IsBand() const { return (attribute & Band) != 0; }
  176.    FREE_CHECK(MatrixType)
  177. };
  178.  
  179. void TestTypeAdd();                            // test +
  180. void TestTypeMult();                           // test *
  181. void TestTypeSP();                             // test SP
  182. void TestTypeOrder();                          // test >=
  183.  
  184.  
  185. /************************* class MatrixBandWidth ***********************/
  186.  
  187. class MatrixBandWidth
  188. {
  189. public:
  190.    int lower;
  191.    int upper;
  192.    MatrixBandWidth(const int l, const int u) : lower(l), upper (u) {}
  193.    MatrixBandWidth(const int i) : lower(i), upper(i) {}
  194.    MatrixBandWidth operator+(const MatrixBandWidth&) const;
  195.    MatrixBandWidth operator*(const MatrixBandWidth&) const;
  196.    MatrixBandWidth minimum(const MatrixBandWidth&) const;
  197.    MatrixBandWidth t() const { return MatrixBandWidth(upper,lower); }
  198.    Boolean operator==(const MatrixBandWidth& bw) const
  199.       { return (lower == bw.lower) && (upper == bw.upper); }
  200.    FREE_CHECK(MatrixBandWidth)
  201. };
  202.  
  203.  
  204. /*********************** Array length specifier ************************/
  205.  
  206. // This class is introduced to avoid constructors such as
  207. //   ColumnVector(int)
  208. // being used for conversions
  209.  
  210. class ArrayLengthSpecifier
  211. {
  212.    int value;
  213. public:
  214.    int Value() const { return value; }
  215.    ArrayLengthSpecifier(int l) : value(l) {}
  216.    FREE_CHECK(ArrayLengthSpecifier)
  217. };
  218.  
  219. /*************************** Matrix routines ***************************/
  220.  
  221.  
  222. class MatrixRowCol;                             // defined later
  223. class MatrixRow;
  224. class MatrixCol;
  225.  
  226. class GeneralMatrix;                            // defined later
  227. class AddedMatrix;
  228. class MultipliedMatrix;
  229. class SubtractedMatrix;
  230. class SPMatrix;
  231. class ConcatenatedMatrix;
  232. class StackedMatrix;
  233. class SolvedMatrix;
  234. class ShiftedMatrix;
  235. class ScaledMatrix;
  236. class TransposedMatrix;
  237. class NegatedMatrix;
  238. class InvertedMatrix;
  239. class RowedMatrix;
  240. class ColedMatrix;
  241. class DiagedMatrix;
  242. class MatedMatrix;
  243. class GetSubMatrix;
  244. class ConstMatrix;
  245. class ReturnMatrixX;
  246. class Matrix;
  247. class nricMatrix;
  248. class RowVector;
  249. class ColumnVector;
  250. class SymmetricMatrix;
  251. class UpperTriangularMatrix;
  252. class LowerTriangularMatrix;
  253. class DiagonalMatrix;
  254. class CroutMatrix;
  255. class BandMatrix;
  256. class LowerBandMatrix;
  257. class UpperBandMatrix;
  258. class SymmetricBandMatrix;
  259. class LinearEquationSolver;
  260. class GenericMatrix;
  261.  
  262.  
  263.  
  264. static MatrixType MatrixTypeUnSp(MatrixType::US);
  265.                         // AT&T needs this
  266.  
  267. class BaseMatrix : public Janitor               // base of all matrix classes
  268. {
  269. protected:
  270.    virtual int search(const BaseMatrix*) const = 0;
  271.                         // count number of times matrix
  272.                         // is referred to
  273. public:
  274.    virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
  275.                         // evaluate temporary
  276.    // for old version of G++
  277.    //   virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
  278.    //   GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
  279. #ifndef TEMPS_DESTROYED_QUICKLY
  280.    AddedMatrix operator+(const BaseMatrix&) const;    // results of operations
  281.    MultipliedMatrix operator*(const BaseMatrix&) const;
  282.    SubtractedMatrix operator-(const BaseMatrix&) const;
  283.    ConcatenatedMatrix operator|(const BaseMatrix&) const;
  284.    StackedMatrix operator&(const BaseMatrix&) const;
  285.    ShiftedMatrix operator+(Real) const;
  286.    ScaledMatrix operator*(Real) const;
  287.    ScaledMatrix operator/(Real) const;
  288.    ShiftedMatrix operator-(Real) const;
  289.    TransposedMatrix t() const;
  290. //   TransposedMatrix t;
  291.    NegatedMatrix operator-() const;                   // change sign of elements
  292.    InvertedMatrix i() const;
  293. //   InvertedMatrix i;
  294.    RowedMatrix AsRow() const;
  295.    ColedMatrix AsColumn() const;
  296.    DiagedMatrix AsDiagonal() const;
  297.    MatedMatrix AsMatrix(int,int) const;
  298.    GetSubMatrix SubMatrix(int,int,int,int) const;
  299.    GetSubMatrix SymSubMatrix(int,int) const;
  300.    GetSubMatrix Row(int) const;
  301.    GetSubMatrix Rows(int,int) const;
  302.    GetSubMatrix Column(int) const;
  303.    GetSubMatrix Columns(int,int) const;
  304. #else
  305.    AddedMatrix& operator+(const BaseMatrix&) const;    // results of operations
  306.    MultipliedMatrix& operator*(const BaseMatrix&) const;
  307.    SubtractedMatrix& operator-(const BaseMatrix&) const;
  308.    ConcatenatedMatrix& operator|(const BaseMatrix&) const;
  309.    StackedMatrix& operator&(const BaseMatrix&) const;
  310.    ShiftedMatrix& operator+(Real) const;
  311.    ScaledMatrix& operator*(Real) const;
  312.    ScaledMatrix& operator/(Real) const;
  313.    ShiftedMatrix& operator-(Real) const;
  314.    TransposedMatrix& t() const;
  315. //   TransposedMatrix& t;
  316.    NegatedMatrix& operator-() const;                   // change sign of elements
  317.    InvertedMatrix& i() const;
  318. //   InvertedMatrix& i;
  319.    RowedMatrix& AsRow() const;
  320.    ColedMatrix& AsColumn() const;
  321.    DiagedMatrix& AsDiagonal() const;
  322.    MatedMatrix& AsMatrix(int,int) const;
  323.    GetSubMatrix& SubMatrix(int,int,int,int) const;
  324.    GetSubMatrix& SymSubMatrix(int,int) const;
  325.    GetSubMatrix& Row(int) const;
  326.    GetSubMatrix& Rows(int,int) const;
  327.    GetSubMatrix& Column(int) const;
  328.    GetSubMatrix& Columns(int,int) const;
  329. #endif
  330.    Real AsScalar() const;                      // conversion of 1 x 1 matrix
  331.    virtual LogAndSign LogDeterminant() const;
  332.    virtual Real SumSquare() const;
  333.    virtual Real SumAbsoluteValue() const;
  334.    virtual Real Sum() const;
  335.    virtual Real MaximumAbsoluteValue() const;
  336.    virtual Real Trace() const;
  337.    Real Norm1() const;
  338.    Real NormInfinity() const;
  339.    virtual MatrixBandWidth BandWidth() const;  // bandwidths of band matrix
  340.    virtual void CleanUp() {}                     // to clear store
  341. //protected:
  342. //   BaseMatrix() : t(this), i(this) {}
  343.  
  344.    friend class GeneralMatrix;
  345.    friend class Matrix;
  346.    friend class nricMatrix;
  347.    friend class RowVector;
  348.    friend class ColumnVector;
  349.    friend class SymmetricMatrix;
  350.    friend class UpperTriangularMatrix;
  351.    friend class LowerTriangularMatrix;
  352.    friend class DiagonalMatrix;
  353.    friend class CroutMatrix;
  354.    friend class BandMatrix;
  355.    friend class LowerBandMatrix;
  356.    friend class UpperBandMatrix;
  357.    friend class SymmetricBandMatrix;
  358.    friend class AddedMatrix;
  359.    friend class MultipliedMatrix;
  360.    friend class SubtractedMatrix;
  361.    friend class SPMatrix;
  362.    friend class ConcatenatedMatrix;
  363.    friend class StackedMatrix;
  364.    friend class SolvedMatrix;
  365.    friend class ShiftedMatrix;
  366.    friend class ScaledMatrix;
  367.    friend class TransposedMatrix;
  368.    friend class NegatedMatrix;
  369.    friend class InvertedMatrix;
  370.    friend class RowedMatrix;
  371.    friend class ColedMatrix;
  372.    friend class DiagedMatrix;
  373.    friend class MatedMatrix;
  374.    friend class GetSubMatrix;
  375.    friend class ConstMatrix;
  376.    friend class ReturnMatrixX;
  377.    friend class LinearEquationSolver;
  378.     friend class GenericMatrix;
  379.    NEW_DELETE(BaseMatrix)
  380. };
  381.  
  382.  
  383. /******************************* working classes **************************/
  384.  
  385. class GeneralMatrix : public BaseMatrix         // declarable matrix types
  386. {
  387.    virtual GeneralMatrix* Image() const;        // copy of matrix
  388. protected:
  389.    int tag;                                     // shows whether can reuse
  390.    int nrows, ncols;                            // dimensions
  391.    int storage;                                 // total store required
  392.    Real* store;                                 // point to store (0=not set)
  393.    GeneralMatrix();                             // initialise with no store
  394.    GeneralMatrix(ArrayLengthSpecifier);         // constructor getting store
  395.    void Add(GeneralMatrix*, Real);              // sum of GM and Real
  396.    void Add(Real);                              // add Real to this
  397.    void Multiply(GeneralMatrix*, Real);         // product of GM and Real
  398.    void Multiply(Real);                         // multiply this by Real
  399.    void Negate(GeneralMatrix*);                 // change sign
  400.    void Negate();                               // change sign
  401.    void operator=(Real);                        // set matrix to constant
  402.    Real* GetStore();                            // get store or copy
  403.    GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
  404.                                                 // temporarily access store
  405.    void GetMatrix(const GeneralMatrix*);        // used by = and initialise
  406.    void Eq(const BaseMatrix&, MatrixType);      // used by =
  407.    int search(const BaseMatrix*) const;
  408.    virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  409.    void CheckConversion(const BaseMatrix&);     // check conversion OK
  410.    void ReDimension(int, int, int);             // change dimensions
  411. public:
  412.    GeneralMatrix* Evaluate(MatrixType);
  413.    virtual MatrixType Type() const = 0;       // type of a matrix
  414.    int Nrows() const { return nrows; }          // get dimensions
  415.    int Ncols() const { return ncols; }
  416.    int Storage() const { return storage; }
  417.    Real* Store() const { return store; }
  418.    virtual ~GeneralMatrix();                    // delete store if set
  419.    void tDelete();                              // delete if tag permits
  420.    Boolean reuse();                                // TRUE if tag allows reuse
  421.    void Protect() { tag=-1; }                   // can't delete or reuse
  422.    int Tag() const { return tag; }
  423.    Boolean IsZero() const;                         // test matrix has all zeros
  424.    void Release() { tag=1; }                    // del store after next use
  425.    void Release(int t) { tag=t; }               // del store after t accesses
  426.    void ReleaseAndDelete() { tag=0; }           // delete matrix after use
  427.    void operator<<(const Real*);                // assignment from an array
  428.    void operator<<(const BaseMatrix& X) { Eq(X,this->Type()); }
  429.                                                 // = without checking type
  430.    void Inject(const GeneralMatrix&);           // copy stored els only
  431.    virtual GeneralMatrix* MakeSolver();         // for solving
  432.    virtual void Solver(MatrixRowCol&, const MatrixRowCol&) {}
  433.    virtual void GetRow(MatrixRowCol&) = 0;      // Get matrix row
  434.    virtual void RestoreRow(MatrixRowCol&) {}    // Restore matrix row
  435.    virtual void NextRow(MatrixRowCol&);         // Go to next row
  436.    virtual void GetCol(MatrixRowCol&) = 0;      // Get matrix col
  437.    virtual void RestoreCol(MatrixRowCol&) {}    // Restore matrix col
  438.    virtual void NextCol(MatrixRowCol&);         // Go to next col
  439.    Real SumSquare() const;
  440.    Real SumAbsoluteValue() const;
  441.    Real Sum() const;
  442.    Real MaximumAbsoluteValue() const;
  443.    LogAndSign LogDeterminant() const;
  444. #ifndef TEMPS_DESTROYED_QUICKLY
  445.    ConstMatrix c() const;                       // to access constant matrices
  446. #else
  447.    ConstMatrix& c() const;                      // to access constant matrices
  448. #endif
  449.    void CheckStore() const;                     // check store is non-zero
  450.    virtual void SetParameters(const GeneralMatrix*) {}
  451.                                                 // set parameters in GetMatrix
  452.    operator ReturnMatrix() const;               // for building a ReturnMatrix
  453.    ReturnMatrix ForReturn() const;
  454.    MatrixInput operator<<(Real);                // for loading a list
  455.    void CleanUp();                              // to clear store
  456.  
  457.    friend class Matrix;
  458.    friend class nricMatrix;
  459.    friend class SymmetricMatrix;
  460.    friend class UpperTriangularMatrix;
  461.    friend class LowerTriangularMatrix;
  462.    friend class DiagonalMatrix;
  463.    friend class CroutMatrix;
  464.    friend class RowVector;
  465.    friend class ColumnVector;
  466.    friend class BandMatrix;
  467.    friend class LowerBandMatrix;
  468.    friend class UpperBandMatrix;
  469.    friend class SymmetricBandMatrix;
  470.    friend class BaseMatrix;
  471.    friend class AddedMatrix;
  472.    friend class MultipliedMatrix;
  473.    friend class SubtractedMatrix;
  474.    friend class SPMatrix;
  475.    friend class ConcatenatedMatrix;
  476.    friend class StackedMatrix;
  477.    friend class SolvedMatrix;
  478.    friend class ShiftedMatrix;
  479.    friend class ScaledMatrix;
  480.    friend class TransposedMatrix;
  481.    friend class NegatedMatrix;
  482.    friend class InvertedMatrix;
  483.    friend class RowedMatrix;
  484.    friend class ColedMatrix;
  485.    friend class DiagedMatrix;
  486.    friend class MatedMatrix;
  487.    friend class GetSubMatrix;
  488.    friend class ConstMatrix;
  489.    friend class ReturnMatrixX;
  490.    friend class LinearEquationSolver;
  491.    friend class GenericMatrix;
  492.    NEW_DELETE(GeneralMatrix)
  493. };
  494.  
  495. class Matrix : public GeneralMatrix             // usual rectangular matrix
  496. {
  497.    GeneralMatrix* Image() const;                // copy of matrix
  498. public:
  499.    Matrix() {}
  500.    Matrix(int, int);                            // standard declaration
  501.    Matrix(const BaseMatrix&);                   // evaluate BaseMatrix
  502.    void operator=(const BaseMatrix&);
  503.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  504.    void operator=(const Matrix& m) { operator=((const BaseMatrix&)m); }
  505.    MatrixType Type() const;
  506.    Real& operator()(int, int);                  // access element
  507.    Real& element(int, int);                     // access element
  508. #ifdef SETUP_C_SUBSCRIPTS
  509.    Real* operator[](int m) { return store+m*ncols; }
  510.    const Real* operator[](int m) const { return store+m*ncols; }
  511. #endif
  512.    Matrix(const Matrix& gm) { GetMatrix(&gm); }
  513. #ifndef __ZTC__
  514.    Real operator()(int, int) const;             // access element
  515. #endif
  516.    GeneralMatrix* MakeSolver();
  517.    Real Trace() const;
  518.    void GetRow(MatrixRowCol&);
  519.    void GetCol(MatrixRowCol&);
  520.    void RestoreCol(MatrixRowCol&);
  521.    void NextRow(MatrixRowCol&);
  522.    void NextCol(MatrixRowCol&);
  523.    virtual void ReDimension(int,int);           // change dimensions
  524.       // virtual so we will catch it being used in a vector called as a matrix
  525.    NEW_DELETE(Matrix)
  526. };
  527.  
  528. class nricMatrix : public Matrix                // for use with Numerical
  529.                                                 // Recipes in C
  530. {
  531.    GeneralMatrix* Image() const;                // copy of matrix
  532.    Real** row_pointer;                          // points to rows
  533.    void MakeRowPointer();                       // build rowpointer
  534.    void DeleteRowPointer();
  535. public:
  536.    nricMatrix() {}
  537.    nricMatrix(int m, int n)                     // standard declaration
  538.       :  Matrix(m,n) { MakeRowPointer(); }
  539.    nricMatrix(const BaseMatrix& bm)             // evaluate BaseMatrix
  540.       :  Matrix(bm) { MakeRowPointer(); }
  541.    void operator=(const BaseMatrix& bm)
  542.       { DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer(); }
  543.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  544.    void operator=(const nricMatrix& m) { operator=((const BaseMatrix&)m); }
  545.    void operator<<(const BaseMatrix& X)
  546.       { DeleteRowPointer(); Eq(X,this->Type()); MakeRowPointer(); }
  547.    nricMatrix(const nricMatrix& gm) { GetMatrix(&gm); MakeRowPointer(); }
  548.    void ReDimension(int m, int n)               // change dimensions
  549.       { DeleteRowPointer(); Matrix::ReDimension(m,n); MakeRowPointer(); }
  550.    ~nricMatrix() { DeleteRowPointer(); }
  551. #ifndef __ZTC__
  552.    Real** nric() const { CheckStore(); return row_pointer-1; }
  553. #endif
  554.    void CleanUp();                                // to clear store
  555.    NEW_DELETE(nricMatrix)
  556. };
  557.  
  558. class SymmetricMatrix : public GeneralMatrix
  559. {
  560.    GeneralMatrix* Image() const;                // copy of matrix
  561. public:
  562.    SymmetricMatrix() {}
  563.    SymmetricMatrix(ArrayLengthSpecifier);
  564.    SymmetricMatrix(const BaseMatrix&);
  565.    void operator=(const BaseMatrix&);
  566.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  567.    void operator=(const SymmetricMatrix& m) { operator=((const BaseMatrix&)m); }
  568.    Real& operator()(int, int);                  // access element
  569.    Real& element(int, int);                     // access element
  570. #ifdef SETUP_C_SUBSCRIPTS
  571.    Real* operator[](int m) { return store+(m*(m+1))/2; }
  572.    const Real* operator[](int m) const { return store+(m*(m+1))/2; }
  573. #endif
  574.    MatrixType Type() const;
  575.    SymmetricMatrix(const SymmetricMatrix& gm) { GetMatrix(&gm); }
  576. #ifndef __ZTC__
  577.    Real operator()(int, int) const;             // access element
  578. #endif
  579.    Real SumSquare() const;
  580.    Real SumAbsoluteValue() const;
  581.    Real Sum() const;
  582.    Real Trace() const;
  583.    void GetRow(MatrixRowCol&);
  584.    void GetCol(MatrixRowCol&);
  585.    void RestoreCol(MatrixRowCol&);
  586.     GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  587.    void ReDimension(int);                       // change dimensions
  588.    NEW_DELETE(SymmetricMatrix)
  589. };
  590.  
  591. class UpperTriangularMatrix : public GeneralMatrix
  592. {
  593.    GeneralMatrix* Image() const;                // copy of matrix
  594. public:
  595.    UpperTriangularMatrix() {}
  596.    UpperTriangularMatrix(ArrayLengthSpecifier);
  597.    void operator=(const BaseMatrix&);
  598.    void operator=(const UpperTriangularMatrix& m)
  599.       { operator=((const BaseMatrix&)m); }
  600.    UpperTriangularMatrix(const BaseMatrix&);
  601.    UpperTriangularMatrix(const UpperTriangularMatrix& gm) { GetMatrix(&gm); }
  602. #ifndef __ZTC__
  603.    Real operator()(int, int) const;             // access element
  604. #endif
  605.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  606.    Real& operator()(int, int);                  // access element
  607.    Real& element(int, int);                     // access element
  608. #ifdef SETUP_C_SUBSCRIPTS
  609.    Real* operator[](int m) { return store+m*ncols-(m*(m+1))/2; }
  610.    const Real* operator[](int m) const { return store+m*ncols-(m*(m+1))/2; }
  611. #endif
  612.    MatrixType Type() const;
  613.    GeneralMatrix* MakeSolver() { return this; } // for solving
  614.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  615.    LogAndSign LogDeterminant() const;
  616.    Real Trace() const;
  617.    void GetRow(MatrixRowCol&);
  618.    void GetCol(MatrixRowCol&);
  619.    void RestoreCol(MatrixRowCol&);
  620.    void NextRow(MatrixRowCol&);
  621.    void ReDimension(int);                       // change dimensions
  622.    NEW_DELETE(UpperTriangularMatrix)
  623. };
  624.  
  625. class LowerTriangularMatrix : public GeneralMatrix
  626. {
  627.    GeneralMatrix* Image() const;                // copy of matrix
  628. public:
  629.    LowerTriangularMatrix() {}
  630.    LowerTriangularMatrix(ArrayLengthSpecifier);
  631.    LowerTriangularMatrix(const LowerTriangularMatrix& gm) { GetMatrix(&gm); }
  632. #ifndef __ZTC__
  633.    Real operator()(int, int) const;             // access element
  634. #endif
  635.    LowerTriangularMatrix(const BaseMatrix& M);
  636.    void operator=(const BaseMatrix&);
  637.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  638.    void operator=(const LowerTriangularMatrix& m)
  639.       { operator=((const BaseMatrix&)m); }
  640.    Real& operator()(int, int);                  // access element
  641.    Real& element(int, int);                     // access element
  642. #ifdef SETUP_C_SUBSCRIPTS
  643.    Real* operator[](int m) { return store+(m*(m+1))/2; }
  644.    const Real* operator[](int m) const { return store+(m*(m+1))/2; }
  645. #endif
  646.    MatrixType Type() const;
  647.    GeneralMatrix* MakeSolver() { return this; } // for solving
  648.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  649.    LogAndSign LogDeterminant() const;
  650.    Real Trace() const;
  651.    void GetRow(MatrixRowCol&);
  652.    void GetCol(MatrixRowCol&);
  653.    void RestoreCol(MatrixRowCol&);
  654.    void NextRow(MatrixRowCol&);
  655.    void ReDimension(int);                       // change dimensions
  656.    NEW_DELETE(LowerTriangularMatrix)
  657. };
  658.  
  659. class DiagonalMatrix : public GeneralMatrix
  660. {
  661.    GeneralMatrix* Image() const;                // copy of matrix
  662. public:
  663.    DiagonalMatrix() {}
  664.    DiagonalMatrix(ArrayLengthSpecifier);
  665.    DiagonalMatrix(const BaseMatrix&);
  666.    DiagonalMatrix(const DiagonalMatrix& gm) { GetMatrix(&gm); }
  667. #ifndef __ZTC__
  668.    Real operator()(int, int) const;             // access element
  669.    Real operator()(int) const;
  670. #endif
  671.    void operator=(const BaseMatrix&);
  672.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  673.    void operator=(const DiagonalMatrix& m) { operator=((const BaseMatrix&)m); }
  674.    Real& operator()(int, int);                  // access element
  675.    Real& operator()(int);                       // access element
  676.    Real& element(int, int);                     // access element
  677.    Real& element(int);                          // access element
  678. #ifdef SETUP_C_SUBSCRIPTS
  679.    Real& operator[](int m) { return store[m]; }
  680.    const Real& operator[](int m) const { return store[m]; }
  681. #endif
  682.    MatrixType Type() const;
  683.  
  684.    LogAndSign LogDeterminant() const;
  685.    Real Trace() const;
  686.    void GetRow(MatrixRowCol&);
  687.    void GetCol(MatrixRowCol&);
  688.    void NextRow(MatrixRowCol&);
  689.    void NextCol(MatrixRowCol&);
  690.    GeneralMatrix* MakeSolver() { return this; } // for solving
  691.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  692.    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  693.    void ReDimension(int);                       // change dimensions
  694. #ifndef __ZTC__
  695.    Real* nric() const
  696.       { CheckStore(); return store-1; }         // for use by NRIC
  697. #endif
  698.    MatrixBandWidth BandWidth() const;
  699.    NEW_DELETE(DiagonalMatrix)
  700. };
  701.  
  702. class RowVector : public Matrix
  703. {
  704.    GeneralMatrix* Image() const;                // copy of matrix
  705. public:
  706.    RowVector() {}
  707.    RowVector(ArrayLengthSpecifier n) : Matrix(1,n.Value()) {}
  708.    RowVector(const BaseMatrix&);
  709.    RowVector(const RowVector& gm) { GetMatrix(&gm); }
  710. #ifndef __ZTC__
  711.    Real operator()(int) const;                  // access element
  712. #endif
  713.    void operator=(const BaseMatrix&);
  714.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  715.    void operator=(const RowVector& m) { operator=((const BaseMatrix&)m); }
  716.    Real& operator()(int);                       // access element
  717.    Real& element(int);                          // access element
  718. #ifdef SETUP_C_SUBSCRIPTS
  719.    Real& operator[](int m) { return store[m]; }
  720.    const Real& operator[](int m) const { return store[m]; }
  721. #endif
  722.    MatrixType Type() const;
  723.    void GetCol(MatrixRowCol&);
  724.    void NextCol(MatrixRowCol&);
  725.    void RestoreCol(MatrixRowCol&);
  726.    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  727.    void ReDimension(int);                       // change dimensions
  728.    void ReDimension(int,int);                   // in case access is matrix
  729. #ifndef __ZTC__
  730.    Real* nric() const
  731.       { CheckStore(); return store-1; }         // for use by NRIC
  732. #endif
  733.    void CleanUp();                              // to clear store
  734.    NEW_DELETE(RowVector)
  735. };
  736.  
  737. class ColumnVector : public Matrix
  738. {
  739.    GeneralMatrix* Image() const;                // copy of matrix
  740. public:
  741.    ColumnVector() {}
  742.    ColumnVector(ArrayLengthSpecifier n) : Matrix(n.Value(),1) {}
  743.    ColumnVector(const BaseMatrix&);
  744.    ColumnVector(const ColumnVector& gm) { GetMatrix(&gm); }
  745. #ifndef __ZTC__
  746.    Real operator()(int) const;                  // access element
  747. #endif
  748.    void operator=(const BaseMatrix&);
  749.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  750.    void operator=(const ColumnVector& m) { operator=((const BaseMatrix&)m); }
  751.    Real& operator()(int);                       // access element
  752.    Real& element(int);                          // access element
  753. #ifdef SETUP_C_SUBSCRIPTS
  754.    Real& operator[](int m) { return store[m]; }
  755.    const Real& operator[](int m) const { return store[m]; }
  756. #endif
  757.    MatrixType Type() const;
  758.    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  759.    void ReDimension(int);                       // change dimensions
  760.    void ReDimension(int,int);                   // in case access is matrix
  761. #ifndef __ZTC__
  762.    Real* nric() const
  763.       { CheckStore(); return store-1; }         // for use by NRIC
  764. #endif
  765.    void CleanUp();                              // to clear store
  766.    NEW_DELETE(ColumnVector)
  767. };
  768.  
  769. class CroutMatrix : public GeneralMatrix        // for LU decomposition
  770. {
  771.    int* indx;
  772.    Boolean d;
  773.    Boolean sing;
  774.    void ludcmp();
  775. public:
  776.    CroutMatrix(const BaseMatrix&);
  777.    MatrixType Type() const;
  778.    void lubksb(Real*, int=0);
  779.    ~CroutMatrix();
  780.    GeneralMatrix* MakeSolver() { return this; } // for solving
  781.    LogAndSign LogDeterminant() const;
  782.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  783.    void GetRow(MatrixRowCol&);
  784.    void GetCol(MatrixRowCol&);
  785.    void operator=(const BaseMatrix&);
  786.    void operator=(const CroutMatrix& m) { operator=((const BaseMatrix&)m); }
  787.    void CleanUp();                                // to clear store
  788.    NEW_DELETE(CroutMatrix)
  789. };
  790.  
  791. /******************************* band matrices ***************************/
  792.  
  793. class BandMatrix : public GeneralMatrix         // band matrix
  794. {
  795.    GeneralMatrix* Image() const;                // copy of matrix
  796. protected:
  797.    void CornerClear() const;                    // set unused elements to zero
  798. public:
  799.    int lower, upper;                            // band widths
  800.    BandMatrix() { lower=0; upper=0; CornerClear(); }
  801.    BandMatrix(int n,int lb,int ub) { ReDimension(n,lb,ub); CornerClear(); }
  802.                                                 // standard declaration
  803.    BandMatrix(const BaseMatrix&);               // evaluate BaseMatrix
  804.    void operator=(const BaseMatrix&);
  805.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  806.    void operator=(const BandMatrix& m) { operator=((const BaseMatrix&)m); }
  807.    MatrixType Type() const;
  808.    Real& operator()(int, int);                  // access element
  809.    Real& element(int, int);                     // access element
  810. #ifdef SETUP_C_SUBSCRIPTS
  811.    Real* operator[](int m) { return store+(upper+lower)*m+lower; }
  812.    const Real* operator[](int m) const { return store+(upper+lower)*m+lower; }
  813. #endif
  814.    BandMatrix(const BandMatrix& gm) { GetMatrix(&gm); }
  815. #ifndef __ZTC__
  816.    Real operator()(int, int) const;             // access element
  817. #endif
  818.    LogAndSign LogDeterminant() const;
  819.    GeneralMatrix* MakeSolver();
  820.    Real Trace() const;
  821.    Real SumSquare() const { CornerClear(); return GeneralMatrix::SumSquare(); }
  822.    Real SumAbsoluteValue() const
  823.       { CornerClear(); return GeneralMatrix::SumAbsoluteValue(); }
  824.    Real Sum() const
  825.       { CornerClear(); return GeneralMatrix::Sum(); }
  826.    Real MaximumAbsoluteValue() const
  827.       { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
  828.    void GetRow(MatrixRowCol&);
  829.    void GetCol(MatrixRowCol&);
  830.    void RestoreCol(MatrixRowCol&);
  831.    void NextRow(MatrixRowCol&);
  832.    void ReDimension(int, int, int);             // change dimensions
  833.    MatrixBandWidth BandWidth() const;
  834.    void SetParameters(const GeneralMatrix*);
  835.    MatrixInput operator<<(Real);                // will give error
  836.    void operator<<(const Real* r);              // will give error
  837.       // the next is included because Zortech and Borland
  838.       // can't find the copy in GeneralMatrix
  839.    void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
  840.    NEW_DELETE(BandMatrix)
  841. };
  842.  
  843. class UpperBandMatrix : public BandMatrix       // upper band matrix
  844. {
  845.    GeneralMatrix* Image() const;                // copy of matrix
  846. public:
  847.    UpperBandMatrix() {}
  848.    UpperBandMatrix(int n, int ubw)              // standard declaration
  849.       : BandMatrix(n, 0, ubw) {}
  850.    UpperBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
  851.    void operator=(const BaseMatrix&);
  852.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  853.    void operator=(const UpperBandMatrix& m)
  854.       { operator=((const BaseMatrix&)m); }
  855.    MatrixType Type() const;
  856.    UpperBandMatrix(const UpperBandMatrix& gm) { GetMatrix(&gm); }
  857.    GeneralMatrix* MakeSolver() { return this; }
  858.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  859.    LogAndSign LogDeterminant() const;
  860.    void ReDimension(int n,int ubw)              // change dimensions
  861.       { BandMatrix::ReDimension(n,0,ubw); }
  862.    Real& operator()(int, int);
  863.    Real& element(int, int);
  864. #ifdef SETUP_C_SUBSCRIPTS
  865.    Real* operator[](int m) { return store+upper*m; }
  866.    const Real* operator[](int m) const { return store+upper*m; }
  867. #endif
  868.    NEW_DELETE(UpperBandMatrix)
  869. };
  870.  
  871. class LowerBandMatrix : public BandMatrix       // upper band matrix
  872. {
  873.    GeneralMatrix* Image() const;                // copy of matrix
  874. public:
  875.    LowerBandMatrix() {}
  876.    LowerBandMatrix(int n, int lbw)              // standard declaration
  877.       : BandMatrix(n, lbw, 0) {}
  878.    LowerBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
  879.    void operator=(const BaseMatrix&);
  880.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  881.    void operator=(const LowerBandMatrix& m)
  882.       { operator=((const BaseMatrix&)m); }
  883.    MatrixType Type() const;
  884.    LowerBandMatrix(const LowerBandMatrix& gm) { GetMatrix(&gm); }
  885.    GeneralMatrix* MakeSolver() { return this; }
  886.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  887.    LogAndSign LogDeterminant() const;
  888.    void ReDimension(int n,int lbw)             // change dimensions
  889.       { BandMatrix::ReDimension(n,lbw,0); }
  890.    Real& operator()(int, int);
  891.    Real& element(int, int);
  892. #ifdef SETUP_C_SUBSCRIPTS
  893.    Real* operator[](int m) { return store+lower*(m+1); }
  894.    const Real* operator[](int m) const { return store+lower*(m+1); }
  895. #endif
  896.    NEW_DELETE(LowerBandMatrix)
  897. };
  898.  
  899. class SymmetricBandMatrix : public GeneralMatrix
  900. {
  901.    GeneralMatrix* Image() const;                // copy of matrix
  902.    void CornerClear() const;                    // set unused elements to zero
  903. public:
  904.    int lower;                                   // lower band width
  905.    SymmetricBandMatrix() { lower=0; CornerClear(); }
  906.    SymmetricBandMatrix(int n, int lb) { ReDimension(n,lb); CornerClear(); }
  907.    SymmetricBandMatrix(const BaseMatrix&);
  908.    void operator=(const BaseMatrix&);
  909.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  910.    void operator=(const SymmetricBandMatrix& m)
  911.       { operator=((const BaseMatrix&)m); }
  912.    Real& operator()(int, int);                  // access element
  913.    Real& element(int, int);                     // access element
  914. #ifdef SETUP_C_SUBSCRIPTS
  915.    Real* operator[](int m) { return store+lower*(m+1); }
  916.    const Real* operator[](int m) const { return store+lower*(m+1); }
  917. #endif
  918.    MatrixType Type() const;
  919.    SymmetricBandMatrix(const SymmetricBandMatrix& gm) { GetMatrix(&gm); }
  920. #ifndef __ZTC__
  921.    Real operator()(int, int) const;             // access element
  922. #endif
  923.    GeneralMatrix* MakeSolver();
  924.    Real SumSquare() const;
  925.    Real SumAbsoluteValue() const;
  926.    Real Sum() const;
  927.    Real MaximumAbsoluteValue() const
  928.       { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
  929.    Real Trace() const;
  930.    LogAndSign LogDeterminant() const;
  931.    void GetRow(MatrixRowCol&);
  932.    void GetCol(MatrixRowCol&);
  933.    void RestoreCol(MatrixRowCol&);
  934.    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  935.    void ReDimension(int,int);                       // change dimensions
  936.    MatrixBandWidth BandWidth() const;
  937.    void SetParameters(const GeneralMatrix*);
  938.    NEW_DELETE(SymmetricBandMatrix)
  939. };
  940.  
  941. class BandLUMatrix : public GeneralMatrix
  942. // for LU decomposition of band matrix
  943. {
  944.    int* indx;
  945.    Boolean d;
  946.    Boolean sing;                                   // TRUE if singular
  947.    Real* store2;
  948.    int storage2;
  949.    void ludcmp();
  950.    int m1,m2;                                   // lower and upper
  951. public:
  952.    BandLUMatrix(const BaseMatrix&);
  953.    MatrixType Type() const;
  954.    void lubksb(Real*, int=0);
  955.    ~BandLUMatrix();
  956.    GeneralMatrix* MakeSolver() { return this; } // for solving
  957.    LogAndSign LogDeterminant() const;
  958.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  959.    void GetRow(MatrixRowCol&);
  960.    void GetCol(MatrixRowCol&);
  961.    void operator=(const BaseMatrix&);
  962.    void operator=(const BandLUMatrix& m) { operator=((const BaseMatrix&)m); }
  963.    NEW_DELETE(BandLUMatrix)
  964.    void CleanUp();                                // to clear store
  965. };
  966.  
  967. /**************************** GenericMatrix class ************************/
  968.  
  969. class GenericMatrix : public BaseMatrix
  970. {
  971.    GeneralMatrix* gm;
  972.    int search(const BaseMatrix* bm) const;
  973.    friend class BaseMatrix;
  974. public:
  975.    GenericMatrix() : gm(0) {}
  976.    GenericMatrix(const BaseMatrix& bm)
  977.       { gm = ((BaseMatrix&)bm).Evaluate(); gm = gm->Image(); }
  978.    GenericMatrix(const GenericMatrix& bm)
  979.       { gm = bm.gm->Image(); }
  980.    void operator=(const GenericMatrix&);
  981.    void operator=(const BaseMatrix&);
  982.    ~GenericMatrix() { delete gm; }
  983.    void CleanUp() { delete gm; gm = 0; }
  984.    void Release() { gm->Release(); } 
  985.    GeneralMatrix* Evaluate(MatrixType) { return gm; };
  986.    MatrixBandWidth BandWidth() const;
  987.    NEW_DELETE(GenericMatrix)
  988. };
  989.  
  990. /***************************** temporary classes *************************/
  991.  
  992. class MultipliedMatrix : public BaseMatrix
  993. {
  994. protected:
  995.    // if these union statements cause problems, simply remove them
  996.    // and declare the items individually
  997.    union { const BaseMatrix* bm1; GeneralMatrix* gm1; };
  998.                           // pointers to summands
  999.    union { const BaseMatrix* bm2; GeneralMatrix* gm2; };
  1000.    MultipliedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  1001.       : bm1(bm1x),bm2(bm2x) {}
  1002.    int search(const BaseMatrix*) const;
  1003.    friend class BaseMatrix;
  1004. public:
  1005.    GeneralMatrix* Evaluate(MatrixType);
  1006.    MatrixBandWidth BandWidth() const;
  1007.    NEW_DELETE(MultipliedMatrix)
  1008. };
  1009.  
  1010. class AddedMatrix : public MultipliedMatrix
  1011. {
  1012. protected:
  1013.    AddedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  1014.       : MultipliedMatrix(bm1x,bm2x) {}
  1015.  
  1016.    friend class BaseMatrix;
  1017. public:
  1018.    GeneralMatrix* Evaluate(MatrixType);
  1019.    MatrixBandWidth BandWidth() const;
  1020. #ifdef __GNUG__
  1021.    void SelectVersion(MatrixType, int&, int&) const;
  1022. #else
  1023.    void SelectVersion(MatrixType, Boolean&, Boolean&) const;
  1024. #endif
  1025.    NEW_DELETE(AddedMatrix)
  1026. };
  1027.  
  1028. class SPMatrix : public AddedMatrix
  1029. {
  1030. protected:
  1031.    SPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  1032.       : AddedMatrix(bm1x,bm2x) {}
  1033.  
  1034.    friend class BaseMatrix;
  1035. public:
  1036.    GeneralMatrix* Evaluate(MatrixType);
  1037.    MatrixBandWidth BandWidth() const;
  1038. #ifdef __GNUG__
  1039.    void SelectVersion(MatrixType, int&, int&) const;
  1040. #else
  1041.    void SelectVersion(MatrixType, Boolean&, Boolean&) const;
  1042. #endif
  1043.    NEW_DELETE(AddedMatrix)
  1044.  
  1045. #ifndef TEMPS_DESTROYED_QUICKLY
  1046.    friend SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
  1047. #else
  1048.    friend SPMatrix& SP(const BaseMatrix&, const BaseMatrix&);
  1049. #endif
  1050. };
  1051.  
  1052. class ConcatenatedMatrix : public MultipliedMatrix
  1053. {
  1054. protected:
  1055.    ConcatenatedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  1056.       : MultipliedMatrix(bm1x,bm2x) {}
  1057.  
  1058.    friend class BaseMatrix;
  1059. public:
  1060.    MatrixBandWidth BandWidth() const;
  1061.    GeneralMatrix* Evaluate(MatrixType);
  1062.    NEW_DELETE(ConcatenatedMatrix)
  1063. };
  1064.  
  1065. class StackedMatrix : public ConcatenatedMatrix
  1066. {
  1067. protected:
  1068.    StackedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  1069.       : ConcatenatedMatrix(bm1x,bm2x) {}
  1070.  
  1071.    friend class BaseMatrix;
  1072. public:
  1073.    GeneralMatrix* Evaluate(MatrixType);
  1074.    NEW_DELETE(StackedMatrix)
  1075. };
  1076.  
  1077. class SolvedMatrix : public MultipliedMatrix
  1078. {
  1079.    SolvedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  1080.       : MultipliedMatrix(bm1x,bm2x) {}
  1081.    friend class BaseMatrix;
  1082.    friend class InvertedMatrix;                        // for operator*
  1083. public:
  1084.    GeneralMatrix* Evaluate(MatrixType);
  1085.    MatrixBandWidth BandWidth() const;
  1086.    NEW_DELETE(SolvedMatrix)
  1087. };
  1088.  
  1089. class SubtractedMatrix : public AddedMatrix
  1090. {
  1091.    SubtractedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  1092.       : AddedMatrix(bm1x,bm2x) {}
  1093.    friend class BaseMatrix;
  1094. public:
  1095.    GeneralMatrix* Evaluate(MatrixType);
  1096.    NEW_DELETE(SubtractedMatrix)
  1097. };
  1098.  
  1099. class ShiftedMatrix : public BaseMatrix
  1100. {
  1101. protected:
  1102.    Real f;
  1103.    union { const BaseMatrix* bm; GeneralMatrix* gm; };
  1104.    ShiftedMatrix(const BaseMatrix* bmx, Real fx) : bm(bmx),f(fx) {}
  1105.    int search(const BaseMatrix*) const;
  1106.    friend class BaseMatrix;
  1107. public:
  1108.    GeneralMatrix* Evaluate(MatrixType);
  1109.    NEW_DELETE(ShiftedMatrix)
  1110. };
  1111.  
  1112. class ScaledMatrix : public ShiftedMatrix
  1113. {
  1114.    ScaledMatrix(const BaseMatrix* bmx, Real fx) : ShiftedMatrix(bmx,fx) {}
  1115.    friend class BaseMatrix;
  1116. public:
  1117.    GeneralMatrix* Evaluate(MatrixType);
  1118.    MatrixBandWidth BandWidth() const;
  1119.    NEW_DELETE(ScaledMatrix)
  1120. };
  1121.  
  1122. class NegatedMatrix : public BaseMatrix
  1123. {
  1124. protected:
  1125.    union { const BaseMatrix* bm; GeneralMatrix* gm; };
  1126.    NegatedMatrix(const BaseMatrix* bmx) : bm(bmx) {}
  1127.    int search(const BaseMatrix*) const;
  1128. private:
  1129.    friend class BaseMatrix;
  1130. public:
  1131.    GeneralMatrix* Evaluate(MatrixType);
  1132.    MatrixBandWidth BandWidth() const;
  1133.    NEW_DELETE(NegatedMatrix)
  1134. };
  1135.  
  1136. class TransposedMatrix : public NegatedMatrix
  1137. {
  1138.    TransposedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1139.    friend class BaseMatrix;
  1140. public:
  1141.    GeneralMatrix* Evaluate(MatrixType);
  1142.    MatrixBandWidth BandWidth() const;
  1143.    NEW_DELETE(TransposedMatrix)
  1144. };
  1145.  
  1146. class InvertedMatrix : public NegatedMatrix
  1147. {
  1148.    InvertedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1149. public:
  1150. #ifndef TEMPS_DESTROYED_QUICKLY
  1151.    SolvedMatrix operator*(const BaseMatrix&) const;  // inverse(A) * B
  1152.     ScaledMatrix operator*(Real t) const { return BaseMatrix::operator*(t); }
  1153. #else
  1154.     SolvedMatrix& operator*(const BaseMatrix&);       // inverse(A) * B
  1155.     ScaledMatrix& operator*(Real t) const { return BaseMatrix::operator*(t); }
  1156. #endif
  1157.    friend class BaseMatrix;
  1158.    GeneralMatrix* Evaluate(MatrixType);
  1159.    MatrixBandWidth BandWidth() const;
  1160.    NEW_DELETE(InvertedMatrix)
  1161. };
  1162.  
  1163. class RowedMatrix : public NegatedMatrix
  1164. {
  1165.    RowedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1166.    friend class BaseMatrix;
  1167. public:
  1168.    GeneralMatrix* Evaluate(MatrixType);
  1169.    MatrixBandWidth BandWidth() const;
  1170.    NEW_DELETE(RowedMatrix)
  1171. };
  1172.  
  1173. class ColedMatrix : public NegatedMatrix
  1174. {
  1175.    ColedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1176.    friend class BaseMatrix;
  1177. public:
  1178.    GeneralMatrix* Evaluate(MatrixType);
  1179.    MatrixBandWidth BandWidth() const;
  1180.    NEW_DELETE(ColedMatrix)
  1181. };
  1182.  
  1183. class DiagedMatrix : public NegatedMatrix
  1184. {
  1185.    DiagedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1186.    friend class BaseMatrix;
  1187. public:
  1188.    GeneralMatrix* Evaluate(MatrixType);
  1189.    MatrixBandWidth BandWidth() const;
  1190.    NEW_DELETE(DiagedMatrix)
  1191. };
  1192.  
  1193. class MatedMatrix : public NegatedMatrix
  1194. {
  1195.    int nr, nc;
  1196.    MatedMatrix(const BaseMatrix* bmx, int nrx, int ncx)
  1197.       : NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
  1198.    friend class BaseMatrix;
  1199. public:
  1200.    GeneralMatrix* Evaluate(MatrixType);
  1201.    MatrixBandWidth BandWidth() const;
  1202.    NEW_DELETE(MatedMatrix)
  1203. };
  1204.  
  1205. class ConstMatrix : public BaseMatrix
  1206. {
  1207.    const GeneralMatrix* cgm;
  1208.    int search(const BaseMatrix*) const;
  1209.    ConstMatrix(const GeneralMatrix* cgmx) : cgm(cgmx) {}
  1210.    friend class BaseMatrix;
  1211.    friend class GeneralMatrix;
  1212. public:
  1213.    GeneralMatrix* Evaluate(MatrixType);
  1214.    MatrixBandWidth BandWidth() const;
  1215.    NEW_DELETE(ConstMatrix)
  1216. };
  1217.  
  1218. class ReturnMatrixX : public BaseMatrix    // for matrix return
  1219. {
  1220.    GeneralMatrix* gm;
  1221.    int search(const BaseMatrix*) const;
  1222. public:
  1223.    GeneralMatrix* Evaluate(MatrixType);
  1224.    friend class BaseMatrix;
  1225. #ifdef TEMPS_DESTROYED_QUICKLY_R
  1226.    ReturnMatrixX(const ReturnMatrixX& tm);
  1227. #else
  1228.    ReturnMatrixX(const ReturnMatrixX& tm) : gm(tm.gm) {}
  1229. #endif
  1230.    ReturnMatrixX(const GeneralMatrix* gmx) : gm((GeneralMatrix*&)gmx) {}
  1231. //   ReturnMatrixX(GeneralMatrix&);
  1232.    MatrixBandWidth BandWidth() const;
  1233.    NEW_DELETE(ReturnMatrixX)
  1234. };
  1235.  
  1236.  
  1237. /**************************** submatrices ******************************/
  1238.  
  1239. class GetSubMatrix : public NegatedMatrix
  1240. {
  1241.    int row_skip;
  1242.    int row_number;
  1243.    int col_skip;
  1244.    int col_number;
  1245.    Boolean IsSym;
  1246.  
  1247.    GetSubMatrix
  1248.       (const BaseMatrix* bmx, int rs, int rn, int cs, int cn, Boolean is)
  1249.       : NegatedMatrix(bmx),
  1250.       row_skip(rs), row_number(rn), col_skip(cs), col_number(cn), IsSym(is) {}
  1251.    GetSubMatrix(const GetSubMatrix& g)
  1252.       : NegatedMatrix(g.bm), row_skip(g.row_skip), row_number(g.row_number),
  1253.       col_skip(g.col_skip), col_number(g.col_number), IsSym(g.IsSym) {}
  1254.    void SetUpLHS();
  1255.    friend class BaseMatrix;
  1256. public:
  1257.    GeneralMatrix* Evaluate(MatrixType);
  1258.    void operator=(const BaseMatrix&);
  1259.    void operator=(const GetSubMatrix& m) { operator=((const BaseMatrix&)m); }
  1260.    void operator<<(const BaseMatrix&);
  1261.    void operator<<(const Real*);                // copy from array
  1262.    void operator=(Real);                        // copy from constant
  1263.    void Inject(const GeneralMatrix&);           // copy stored els only
  1264.    MatrixBandWidth BandWidth() const;
  1265.    NEW_DELETE(GetSubMatrix)
  1266. };
  1267.  
  1268. /**************************** matrix input *******************************/
  1269.  
  1270. class MatrixInput          // for reading a list of values into a matrix
  1271.                            // the difficult part is detecting a mismatch
  1272.                            // in the number of elements
  1273. {
  1274.    static int n;           // number values still to be read
  1275.    static Real* r;         // pointer to last thing read
  1276.    static int depth;       // number of objects of this class in existence
  1277. public:
  1278.    MatrixInput() { depth++; }
  1279.    MatrixInput(const MatrixInput&) { depth++; }
  1280.    ~MatrixInput();
  1281.    MatrixInput operator<<(Real);
  1282.                            // could return a reference if we don't have
  1283.                            // TEMPS_DESTROYED_QUICKLY
  1284.    friend class GeneralMatrix;                             
  1285. };
  1286.  
  1287. /***************************** exceptions ********************************/
  1288.  
  1289. class MatrixDetails
  1290. {
  1291.    MatrixType type;
  1292.    int nrows;
  1293.    int ncols;
  1294.    int ubw;
  1295.    int lbw;
  1296. public:
  1297.    MatrixDetails(const GeneralMatrix& A);
  1298.    void PrintOut();
  1299. };
  1300.    
  1301.  
  1302.  
  1303. class SpaceException : public Exception
  1304. {
  1305. public:
  1306.    static long st_type() { return 2; }
  1307.    long type() const { return 2; }
  1308.    static int action;
  1309.    SpaceException();
  1310.    static void SetAction(int a) { action=a; }
  1311.    SpaceException(const SpaceException&) {}
  1312. };
  1313.  
  1314.  
  1315. class MatrixException : public Exception
  1316. {
  1317. protected:
  1318.    MatrixException() {}
  1319. public:
  1320.    static long st_type() { return 3; }
  1321.    long type() const { return 3; }
  1322.    MatrixException(int);
  1323.    MatrixException(int, const GeneralMatrix&);
  1324.    MatrixException(int, const GeneralMatrix&, const GeneralMatrix&);
  1325.    MatrixException(const MatrixException&) {}
  1326. };
  1327.  
  1328. class DataException : public MatrixException
  1329. {
  1330. protected:
  1331.    DataException() {}
  1332. public:
  1333.    static long st_type() { return 3*53; }
  1334.    long type() const { return 3*53; }
  1335.    static int action;
  1336.    DataException(const GeneralMatrix& A);
  1337.    static void SetAction(int a) { action=a; }
  1338.    DataException(const DataException&) {}
  1339. };
  1340.  
  1341. class SingularException : public DataException
  1342. {
  1343. public:
  1344.    static long st_type() { return 3*53*109; }
  1345.    long type() const { return 3*53*109; }
  1346.    SingularException(const GeneralMatrix& A);
  1347.    SingularException(const SingularException&) {}
  1348. };
  1349.  
  1350. class NPDException : public DataException     // Not positive definite
  1351. {
  1352. public:
  1353.    static long st_type() { return 3*53*113; }
  1354.    long type() const { return 3*53*113; }
  1355.    NPDException(const GeneralMatrix&);
  1356.    NPDException(const NPDException&) {}
  1357. };
  1358.  
  1359. class ConvergenceException : public MatrixException
  1360. {
  1361. public:
  1362.    static long st_type() { return 3*59; }
  1363.    long type() const { return 3*59; }
  1364.    static int action;
  1365.    ConvergenceException(const GeneralMatrix& A);
  1366.    static void SetAction(int a) { action=a; }
  1367.    ConvergenceException(const ConvergenceException&) {}
  1368. };
  1369.  
  1370. class ProgramException : public MatrixException
  1371. {
  1372. protected:
  1373.    ProgramException() {}
  1374.    ProgramException(int);            // int is not used
  1375. public:
  1376.    static long st_type() { return 3*61; }
  1377.    long type() const { return 3*61; }
  1378.    static int action;
  1379.    ProgramException(char* c);
  1380.    ProgramException(char* c, const GeneralMatrix&);
  1381.    ProgramException(char* c, const GeneralMatrix&, const GeneralMatrix&);
  1382.    ProgramException(char* c, MatrixType, MatrixType);
  1383.    ProgramException(const GeneralMatrix&);
  1384.    static void SetAction(int a) { action=a; }
  1385.    ProgramException(const ProgramException&) {}
  1386. };
  1387.  
  1388. class IndexException : public ProgramException
  1389. {
  1390. public:
  1391.    static long st_type() { return 3*61*101; }
  1392.    long type() const { return 3*61*101; }
  1393.    IndexException(int i, const GeneralMatrix& A);
  1394.    IndexException(int i, int j, const GeneralMatrix& A);
  1395.    // next two are for access via element function
  1396.    IndexException(int i, const GeneralMatrix& A, Boolean);
  1397.    IndexException(int i, int j, const GeneralMatrix& A, Boolean);
  1398.    IndexException(const IndexException&) {}
  1399. };
  1400.  
  1401. class VectorException : public ProgramException    // can't convert to vector
  1402. {
  1403. public:
  1404.    static long st_type() { return 3*61*107; }
  1405.    long type() const { return 3*61*107; }
  1406.    VectorException();
  1407.    VectorException(const GeneralMatrix& A);
  1408.    VectorException(const VectorException&) {}
  1409. };
  1410.  
  1411. class NotSquareException : public ProgramException
  1412. {
  1413. public:
  1414.    static long st_type() { return 3*61*109; }
  1415.    long type() const { return 3*61*109; }
  1416.    NotSquareException(const GeneralMatrix& A);
  1417.    NotSquareException(const NotSquareException&) {}
  1418. };
  1419.  
  1420. class SubMatrixDimensionException : public ProgramException
  1421. {
  1422. public:
  1423.    static long st_type() { return 3*61*113; }
  1424.    long type() const { return 3*61*113; }
  1425.    SubMatrixDimensionException();
  1426.    SubMatrixDimensionException(const SubMatrixDimensionException&) {}
  1427. };
  1428.  
  1429. class IncompatibleDimensionsException : public ProgramException
  1430. {
  1431. public:
  1432.    static long st_type() { return 3*61*127; }
  1433.    long type() const { return 3*61*127; }
  1434.    IncompatibleDimensionsException();
  1435.    IncompatibleDimensionsException(const IncompatibleDimensionsException&) {}
  1436. };
  1437.  
  1438. class NotDefinedException : public ProgramException
  1439. {
  1440. public:
  1441.    static long st_type() { return 3*61*131; }
  1442.    long type() const { return 3*61*131; }
  1443.    NotDefinedException(char* op, char* matrix);
  1444.    NotDefinedException(const NotDefinedException&) {}
  1445. };
  1446.  
  1447. class CannotBuildException : public ProgramException
  1448. {
  1449. public:
  1450.    static long st_type() { return 3*61*137; }
  1451.    long type() const { return 3*61*137; }
  1452.    CannotBuildException(char* matrix);
  1453.    CannotBuildException(const CannotBuildException&) {}
  1454. };
  1455.  
  1456.  
  1457. class InternalException : public MatrixException
  1458. {
  1459. public:
  1460.    static long st_type() { return 3*67; }
  1461.    long type() const { return 3*67; }
  1462.    static int action;
  1463.    InternalException(char* c);
  1464.    static void SetAction(int a) { action=a; }
  1465.    InternalException(const InternalException&) {}
  1466. };
  1467.  
  1468.  
  1469. /***************************** functions ***********************************/
  1470.  
  1471.  
  1472. inline LogAndSign LogDeterminant(const BaseMatrix& B)
  1473.    { return B.LogDeterminant(); }
  1474. inline Real SumSquare(const BaseMatrix& B) { return B.SumSquare(); }
  1475. inline Real Trace(const BaseMatrix& B) { return B.Trace(); }
  1476. inline Real SumAbsoluteValue(const BaseMatrix& B)
  1477.    { return B.SumAbsoluteValue(); }
  1478. inline Real Sum(const BaseMatrix& B)
  1479.    { return B.Sum(); }
  1480. inline Real MaximumAbsoluteValue(const BaseMatrix& B)
  1481.    { return B.MaximumAbsoluteValue(); }
  1482. inline Real Norm1(const BaseMatrix& B) { return B.Norm1(); }
  1483. inline Real Norm1(RowVector& RV) { return RV.MaximumAbsoluteValue(); }
  1484. inline Real NormInfinity(const BaseMatrix& B) { return B.NormInfinity(); }
  1485. inline Real NormInfinity(ColumnVector& CV)
  1486.    { return CV.MaximumAbsoluteValue(); } 
  1487.  
  1488. #endif
  1489.